function result = phi(t, r)
%PHI This function gets 2 arguments (t,r) and returns phi(t,r)
%   This function implements the following formula:
%   phi(t ,r) = -C * integral(-inf => -epsilon)[G(t,r, -s') * e^(i*omega*s') ds']
%               +C * integral(epsilon => inf)[G(t,r, -s') * e^(i*omega*s') ds']

% Use global parameters:
global C;
global eps;
global omega;
global inf_s;
global h_s;

DBG('disp([''phi(t='' num2str(t) '', r='' num2str(t)])');
%disp(['phi(t=' num2str(t) ', r=' num2str(r) ')']);

%integrand = inline('G_(t, r, s_) .* exp(1i .* omega .* s_)', 't', 'r', 's_', 'omega');
%integrand = inline('G_(t, r, s_) .* exp(1i .* omega .* s_)', 't', 'r', 's_', 'omega');
%integrand = inline('G_(t, r, -500) .* exp(1i .* omega .* s_)', 't', 'r', 's_', 'omega');

low = -C .* trpzInt(@(s_)int_s(t, r, s_, omega), -inf_s, -eps, h_s);  
high =  C .* trpzInt(@(s_)int_s(t, r, s_, omega), eps, inf_s, h_s);  
result = low + high;

end

